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We propose a model for charge density waves in ring shaped crystals, which depicts frustration 
between intra- and inter-chain couplings coming from cylindrical bending. It is then mapped to a 
three dimensional uniformly frustrated XY model with one dimensional anisotropy in connectivity. 
The nonequilibrium relaxation dynamics is investigated by Monte Carlo simulations to find a phase 
transition which is quite different from that of usual whisker crystal. We also find that the low 
temperature state is a three dimensional phase vortex lattice with a two dimensional phase coherence 
in a cylindrical shell and the system shows power law relaxation in the ordered phase. 
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I. INTRODUCTION 



Transition metal chalcogenides such as NbSe3 and 
TaS3 are known as quasi one dimensional materials. Low- 
ering temperature, they take a phase transition to the 
charge density wave (CDW) phase where both of atom 
positions and electronic charge density are modulated in 
twice of the Fermi wave number pj. Recently, Tanda 
and Tsuneta et al. have succeeded to synthesize various 
types of single crystals in closed loop shapes, e.g., simple 
ring, Mobius ring and figure-of-eight, of these materials 
by the chemical vapor transportation methodQ, Ol- X- 
ray diffraction measurement shows that they have the 
same crystalline structure with usual whisker crystals 
and the temperature dependence of conductivity indi- 
cates a CDW transition at the temperature slightly lower 
than the critical temperature of whisker crystals [j] . The 
influence of the crystal geometry or topology on the prop- 
erty is an interesting problem. 



Recently, Shimatake and Toda investigated nonequi- 
librium relaxation dynamics of NbSe3 by using ultrafast 
laser |5j ■ They measured both of ring and whisker crys- 
tals with comparable dimensions and found a significant 
difference between them in the low temperature CDW 
phase. The relaxation time for the initial rapid decay 
shows divergent behavior at the transition temperature 
for a whisker crystal but does not for a ring crystal, ft 
suggests that the phase transition disappears or the type 
of a phase transition changes for ring crystals. In this ar- 
ticle, we propose a simple model for CDWs in ring crys- 
tals and analyze the nonequilibrium relaxation dynamics 
by Monte Carlo method to understand the reason for 
such difference. 
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FIG. 1: Schematic diagram of ring crystal and mapping onto 
cuboid. 



II. MODEL FOR CDW IN RING CRYSTAL 

At first, we introduce a phase field model for CDW in 
a ring crystal. The crystal axes are named a-, b-, and c- 
axes as illustrated in Figure^ The quasi one dimensional 
chains run along the b-axis. The position in a sample is 
expressed as r = (x a , Xb, x c ) — (r,rcf),x c ) by using cylin- 
drical coordinates. The modulated part of charge density 
can be expressed as p(r, t) — p cos (Q(r)xb + 6(r, t)) in 
each chain. Here, Q(r) is the mean wave number of a 
chain at radius r and 6(r, t) is the phase fluctuation vari- 
able. Periodic boundary condition along a chain yields 
Q(r)2nr = 27r N w (r) where N w (r) is a number of waves 
in a chain. Q(r) possibly deviates from the natural wave 
number of CDW, Qo = 2k p, to synchronize the charge 
density modulations between neighboring chains. Such 
commensuration between the chains neighboring on the 
radial direction results strain. The stress grows with the 
thickness and is released by making some imperfection, 
such as a jump of N w , by every several chains. We con- 
sider that the characteristic relaxation dynamics of a ring 
crystal is due to such conflict between two requirements 
for the natural wave length inside each chain and for pe- 
riod matching between neighboring chains. 

Here we suppose bundles of looped chains in which 
no defect exists and then N w (r) is constant. The mean 
wave number decreases with radius as Q(r) = N w /r in 
each bundle. The wave number is larger than the natu- 
ral one for inner edges and smaller for outer edges. The 
Coulomb interaction V a at the boundary of neighboring 
two bundles is proportional to the product of charge den- 
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sities as 

V a oc J dx b cos[QiXb + 6±(xb)] cos{(Qtx b + 9^(x b )] 

= ^ J dxb{cos[9i(x b ) - 6» T (x 6 ) - AQx b ] 
+ coa[B i (x b ) + O^Xb) + (Q l + Q T )x & ]} (1) 

The down- (J,) and up- (f) arrows indicate quantities of 
the inner and outer bundles, respectively. AQ(r) = Qf — 
Ql is a gap of wave number at the boundary between 
the two neighboring bundles. In this expression we keep 
track of the first term only and neglect the second term 
because it oscillates much faster than 9(xb) changes. 

Knowing the above interaction form between the 
phases 0j and 9±, let us consider a discrete version of 
the phase Hamiltonian. We map the ring crystal onto 
a simple cubic lattice with size N = L a x x L c (See 
Figure ^) . Note that the space is divided uniformly in 
cj> instead of x b — re/). The lattice spacings are taken 
as a mesoscopic scale, much larger than the CDW wave 
length. Adding intra-chain distortion energy and inter- 
bundle energy on the c-axis to eq. , the Hamiltonian is 
written as 



H 



'i+a 



+&) 



(2) 



Here, i labels the lattice points of a simple cubic lattice 
and a is the unit lattice vector on the a-axis. What is 
important is that there is an additional term, 



Ay+a = AQx b 5 aa = AN U 



2nfi b 6 aa , (3) 



which comes from the ring geometry. AN W — rAQ is a 
gap of a number of waves and / = AN W /L ( j 1 is a filling 
factor that denotes the mean density of vortices. Ay 
causes a frustration between the couplings along a- and 
b-axes and yields a constant number of phase vortex lines 
directed to the c-axis even without thermal excitation. 
The vorticity is given by the loop integration of V# along 
each plaquette. By omitting Ay, we obtain a usual XY 
model, which is applicable to CDWs in whisker crystals. 

The Hamiltonian eq.J2J) is equivalent to the uniformly 
frustrated XY model which is proposed for the super- 
conductors under magnetic field parallel to the c-axis 
0, 0, 0, where Ay is translated to a vector potential 
in Landau gauge. The anisotropy is, however, quite dif- 
ferent. While high T c superconductor has strong cou- 
pling inside each CuC>2 plane and weak one inter-planes, 
CDW materials has quasi one dimensional property. The 
anisotropy in coupling constants is taken as J b = 7 J a — 
jJ c ,7 > 1 in the present model. Since the distance be- 
tween lattice points neighboring on the b-axis is propor- 
tional to r, coupling constants and filling factor should 
depend on the radius r but we ignore this r-dependence 
as a first approximation. 

At first, we investigate the ground state of this Hamil- 
tonian. By simulated annealing, we obtained an energy 



minimal state for sufficiently large 7. This state has AN W 
phase vortex lines between every neighboring bc-planes, 
which are straight and parallel to the c-axis. These vor- 
tex lines form a two dimensional lattice in the ab-plane 
with the unit lattice vectors (±ao, fro/2/, 0). Order pa- 
rameter r v , which is related to the Bragg peak height of 
a vortex lattice, can be defined as follows; 

r„ = S v (q v )/N, q v = (2vr/2ao, 27r//fr , 0), 
S v (q) = A-^Mq)! 2 ), 

«c(q) = £(vj •c)exp(iq-r j ). (4) 



Here (vj • c) is the vorticity defined at the center of a pla- 
quette j which is perpendicular to the c-axis. Owing to 
the strong coupling along the b-axis, the phase is almost 
uniform along the direction but gently modulated with 
a period bo/ f. This state looks quite different from the 
solitonic solution in the isotropic case where distortion 
is localized around the vortex. Simple variational calcu- 
lation supposing the periods 2ao and fro// along a- and 
b-axes, respectively, yields an energy minimal solution 
expressed as 



% = (-iy 



27- 



■ cos(27r fi b ) 



+ 0(7~ 2 ), (5) 



which agrees well with the simulated annealing result. 
Note that there are many energetically degenerated 
states which are obtained by the transformation: 9i — > 
9\ + 27tni a /L a ,(n — ±1, ±2, •••). This transformation 
does not change the positions of vortices. 



III. NONEQUILIBRIUM RELAXATION 
ANALYSIS OF PHASE TRANSITION 

Next we investigate the relaxation dynamics from an 
ordered state to an equilibrium state at finite temper- 
ature, i.e., a rapid heating process. At the same time, 
the phase transition is also analyzed by the nonequi- 
librium relaxation method [9|- In order to obtain an 
initial state, we set the phases as eq.© and make the 
system relax at zero temperature in diffusive dynamics. 
After that, Metropolis dynamics at finite temperature is 
started. Filling factor / is fixed to 1/16 in this work. 
The same calculation is done for several anisotropy pa- 
rameters, 7 = 10, 16, 24, 32 and 64. The Monte Carlo flip 
is repeated up to 65,000 steps per each site at maximum. 
This step is identified as time, t. Average is taken over 
8-16 samples for each temperature. The sample size used 
is L a x Lc/, x L c = 64 x 512 x 64 and periodic boundary 
condition is imposed for all directions. For this system 
size, finite size effect is not crucial within the observation 
time. 

The order parameter r v equals f 2 in the initial state 
and decreases to the equilibrium value which is zero when 
T > T c and finite when T < T c . To find the critical point, 
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FIG. 2: The local exponent of the order parameter r v for 
various temperatures, as a function of 1/t. The temperature 
of each curve is shown in the figure in the unit of J a /ks (T c — 
6.87J a /fcs)- The derivative is calculated after local fitting to a 
4th order polynomial. The anisotropy parameter 7 equals 32. 
The horizontal line shows A c obtained by dynamical scaling. 

we calculate the local exponent Q, 

X(t) = -d\nr v (t)/d\nt. (6) 

In Figure El A(t) is plotted with respect to 1/t. In the 
limit of t — ► 00, X(t) goes to infinity for T > T c and to 
zero for T < T c . Just at the critical temperature, A(i) 
converges to a certain finite exponent A c which charac- 
terize the critical power law decay. 

For T > T c , r v decays in power function of t for short 
time scale and makes a crossover to exponential decay. 
We can perform dynamical scaling as 

r v (t)=T(T)- x °S v (t/T(T)), (7) 

Here t(T) is a relaxation time, which diverges at the 
critical temperature T c as 

r(T) oc (T/T c - \)~ zv . (8) 

We initially obtain A c and t(T)'s by scaling to eq.Ql then 
fitting of t(T)'s to eq.JSJ yields T c and zv. The scaling 
result for 7 = 32 is shown in Figure^ The results for sev- 
eral 7's are summarized in tabled T c is proportional to 
(J b J c ) 1/2 , which measures the effective elasticity of vor- 
tex lines confined between two bc-planes. The exponent 
A c seems not universal among different 7's. It shows, 
however, a sign for saturation to the anisotropic limit 
value w 0.25 as 7 becomes large, zv is less dependent of 
7 even for small 7. 

For small 7 dynamical scaling does not work well in 
the very vicinity of the critical temperature, where the 
functional form of relaxation changes with temperature. 
In Figure the local exponent for 7 = 10 is plotted 
with 1/t. If the relaxation form is simple exponential 
one in the long time limit, X(t) would be asymptotically 
proportional to t. Figure 0] shows, however, X(t) oc i M 




t/x(T) 

FIG. 3: Scaling plot for 7 = 32. The temperatures for used 
data are 6.9, 7.0, 7.1, 7.2, 7.5, 8.0, 8.5, 9.0, 10.0, 11.0, 12.0 and 
14.0J a /fcs. The inset shows critical behavior of the relaxation 
time. 



where /i is the larger than unity and increases up to 2 as 
approaching the critical temperature. Such bad scaling 
region in temperature, however, becomes narrower as 7 
increases. 



IV. RELAXATION DYNAMICS IN ORDERED 
PHASE 

Next, we investigate the relaxation function in the low 
temperature phase. In order to eliminate the influence 
of the finite equilibrium value, the time derivative of the 
order parameter is calculated. The relaxation of r v does 
not look like an exponential function in the time range 
of the present simulations. We don't find an apparent 
reduction of relaxation time when leaving from the criti- 
cal point unlike in the high temperature phase. Although 
the functional form of r v is unclear, the power law behav- 
ior is observed to be suitable for some other quantities, 
e.g., total energy of the system and the order parameter 
of phase periodicity r h s . Here 

r b e ee S e (4)/N. (9) 



7 


T c (LE) 


T c (DS) 


A (DS) 


ZV (DS) 


10 


1.159(2) 


1.13 


0.06 


2.2 


16 


1.190(5) 


1.17 


0.16 


2.1 


24 


1.208(6) 


1.18 


0.22 


2.2 


32 


1.214(5) 


1.20 


0.24 


2.2 


64 


1.206(6) 


1.18 


0.25 


2.6 



TABLE I: Critical temperature in the unit of ( Jt J c ) 1//2 /fcB 
and exponents for several anisotropy parameters. LE means 
the critical temperature is estimated from the long time be- 
havior of local exponent and DS means dynamical scaling. 
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FIG. 4: The local exponent for a relatively small anisotropic 
parameter 7 = 10. 

Here Sg(q) is a Fourier transformation of phase correla- 
tion function defined as 

5 e (q) = N- 1 c os(0i - 6j) exp [iq • (r ; - r,)] , (10) 

which has peaks at = (2-7r/a ,0, 0) and q# = 
(0, 27r//6o, 0). Note that this order parameter is finite 
at t = only because the initial state is chosen as eq.© 
and other metastable states have peaks in different wave 
number vectors. The time evolution of r b s and tdrg/dt is 
shown in Figure |SJ For T < T c , r b g seems to decay to a 
finite value in the limit of t — > 00. This indicates that not 
only vortex but also phase has true long range periodic 
order in three dimensions. We confirm that this finite 
value does not come from the finite size effect. The time 
derivative shows power law relaxation with t. The decay 
exponent becomes larger with T. 

V. DISCUSSION 
A. dimensionality of periodic order 

So far, it is shown that the quasi one dimensional frus- 
trated XY model takes a vortex lattice melting transition. 
The dynamical scaling analysis strongly suggests that it 
is a standard second order transition. Then r v and r h s 
work as order parameters of vortex and phase itself. 

By equilibrium simulations, we find that hclicity mod- 
ulus, which measures the thermodynamic stiffness of 
phase against the perturbation, — > 6>i + ei a (e — ► 0) , 
is zero along the a-axis for all temperature while those 
along b- and c-axes are finite below T c . This means that 
the (quasi) long range phase coherence is established only 
in the ab-plane. 

The reason why the two dimensional phase coherence 
appears in quasi one dimensional system is considered as 
follows: Even with J a = J C (<C Jb), the a- and c-axes are 




FIG. 5: Time evolution of phase structure factor r\ (top) and 
time derivative of r g multiplied by t (bottom). The derivative 
is obtained by local fitting to a 4th order cubic polynomial. 
The line of slope -0.22 means the relaxation at the critical 
temperature. 



not equivalent since the former conflicts with the b-axis 
and the latter does not. The fast growth of the order 
along the b-axis suppresses the order along the a-axis 
then two dimensional-like order appears. 

Power law behavior and two dimensional order remind 
us of the Kosterlitz-Thouless (KT) phase ^3 with quasi 
long range order. If so, it is strange that r h g is finite in the 
long time limit as shown in the previous section. Further- 
more finite r h e seems directly inconsistent with the fact 
that helicity modulus for the a-axis is zero. Of course 
we can not eliminate the possibility that r h g decreases 
to zero in the time scale much longer than the present 
observation time but there can be another explanation 
as follows. When the two dimensional true long range 
order is established, the center of mass of its phase is 
individually pinned in every bc-plane due to the break 
of ergodicity. As a result the phase fluctuation along 
the a-axis is suppressed even though there is no restoring 
force. Thus finite r v which is owing to the initial condi- 
tion survives even though long range order exists only in 
two dimensions. In contrast, vorticity is a gauge invari- 
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ant quantity and really has three dimensional long range 
order. If this is true, there would be another mechanism 
for the power decay which differs from the KT scenario, 
on which r\ goes to zero. Indeed, power divergence of 
relaxation time at T c is inconsistent with the exponential 
divergence for the KT transition 

B. comparison with in-plane isotropic case 

Here we compare the present result with high T c super- 
conductors. The case of the magnetic field perpendicular 
to the CuC>2 plane is given by J a = J& ^> J c in the present 
notation. This is isotropic in the frustrated ab-plane. It 
is believed that there is a first order transition [j| into 
the low temperature phase where phase coherence exists 
only along the c-axis 7]. The limit 7 — » 1 in the present 
model is also expected to show the same transition al- 
though the discontinuous property is quite small in the 
case that J c is as large as J a , Jb 0- There is a possibility 
that bad scaling behavior near T c for smaller 7 as men- 
tioned at the last of section UTTl is due to the crossover to 
a first order transition. 

Power law relaxation is also observed in the isotropic 
frustration model [6j . This is apparently unrelated to the 
KT order. We consider that the power law behaviors of 
the isotropic and anisotropic models is a common prop- 
erty of the partially ordered states of frustrated systems 
although dimensionality of phase ordering is different. 

C. comparison with quasi two dimensional case 

The case of a magnetic field parallel to the CuC>2 plane 
can be expressed by the present model with 7 <C 1 al- 
though the CuC>2 plane is normal to the b-axis, not to 
the c-axis 13]. This is just the opposite anisotropy to 
the present study and the system is divided into isolated 
planes without frustration in the strong anisotropy limit. 
This case, however, has many common properties with 
the present case of 7 ^ 1 in spite of the difference in 
quasi two and one dimensional properties, e.g., two di- 
mensional phase coherence parallel to the c-axis and the 
same structure of vortex lattices (HIH. Additionally 
the phase transitions change from first order to second 
order as anisotropy becomes stronger in both cases. This 
similarity comes from the fact that quantitative behav- 
ior is not affected by the strength of the coupling along 
the c-axis, which is not concerned with the frustration 
of the system. It is essential that both of these opposite 



anisotropies, 7 ^> 1 and 7 < 1, break the balance of 
frustration in the ab-plane in the same way. 

The vortex lattice order in the quasi two dimensional 
system is considered to be of quasi long range 01 • ^ 
is also reported that the two layer system (this is purely 
two dimensional) has two types of KT phases [Pa. [l5| , 
the vortex phase and the Meissner phase. It is, how- 
ever, not clear whether the inter layer coupling in three 
dimensional system could be irrelevant or not for small 
but finite 7. We consider that true long range order as- 
sisted by three dimensional vortex lattice order is worth 
to re-examine in quasi two dimensional case as well as 
the present case. 

D. Experiment 

Finally let us mention about the relation to the exper- 
iments of a CDW in a ring crystal 5] . The present work 
based on the anisotropic frustrated XY model found a 
phase transition to the ordered state where phase vor- 
tex lines along the c-axis form a lattice and two dimen- 
sional phase coherence is established in each cylindrical 
shell perpendicular to radius direction. Thus it suggests 
that the ring geometry makes the CDW phase transition 
quite different from that of a whisker crystal. The lat- 
ter is described by the XY model without frustration, 
which shows a simple ferromagnetic transition. In the 
low temperature phase of a ring crystal, the relaxation 
of phase fluctuation, which is closely connected to elec- 
tric polarization and reflectivity, shows a power law re- 
laxation without characteristic time scale, so that small 
apparent relaxation time, which does not show singular 
behavior, would be estimated if one supposes an expo- 
nential decay as for whisker crystals. 

There remain some future works. In this work, rapid 
heating process is studied because of its self averaging 
property instead of rapid cooling process. The latter is a 
challenging work because it is proper for the direct com- 
parison with the laser pumping experiment and contains 
the process to untwist entangled flux lines. The present 
model drops some features of the system, e.g., radius de- 
pendence of the model parameters, which causes distri- 
bution of local critical temperature and relaxation time 
then the transition could be blurred, and a pinning effect 
of lattice defects in atomic crystal. 

The authors thank S. Tanda, Y. Toda and K. Shi- 
matake for useful discussions. This work is supported 
by 21st Century COE program "Topological Science and 
Technology" . 



[1] G. Griiner, Rev. Mod. Phys. 60, 1129 (1988). (2004). 

[2] S. Tanda, T. Tsuneta, Y. Okajima, K. Inagaki, K. Ya- [4] S. Tsuneta, S. Tanda, Y. Okajima, K. Inagaki, and K. Ya- 

maya, and N. Hatakenaka, Nature 417, 397 (2002). maya, Physica B 329, 1544 (2003). 

[3] T. Tsuneta and S. Tanda, J. Crystal Growth 264, 223 [5] K. Shimatake, Y. Toda, and S. Tanda, cond- 



6 



mat/0511318 (unpublished). 
[6] Y. H. Li and S. Teitel, Phys. Rev. B 47, 359 (1993). 
[7] T. Chen and S. Teitel, Phys. Rev. B 55, 11766 (1997). 
[8] X. Hu, S. Miyashita, and M. Tachiki, Phys. Rev. Lett. 

79, 3498 (1997). 
[9] Y. Ozcki, K. Ogawa, and N. Ito, Phys. Rev. E 67, 026702 

(2003). 

[10] J. M. Kosterlitz and D. J. Thouless, J. Phys. C 6, 1181 
(1973). 



[11] X. Hu and M. Tachiki, Phys. Rev. Lett. 85, 2577 (2000). 
[12] X. Hu and M. Tachiki, Phys. Rev. B 70, 064506 (2004). 
[13] L. Balents and D. R. Nelson, Phys. Rev. Lett. 73, 2618 
(1994). 

[14] E. Orignac and T. Giamarchi, Phys. Rev. B 64, 144515 
(2001). 

[15] E. Granato, Phys. Rev. B 72, 104521 (2005). 



